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Abstract. Based on a local approximation of the Riemannian distance on a manifold by a computationally 
cheap dissimilarity measure, a time discrete geodesic calculus is developed, and applications to shape space are 
explored. The dissimilarity measure is derived from a deformation energy whose Hessian reproduces the underlying 
Riemannian metric, and it is used to define length and energy of discrete paths in shape space. The notion of discrete 
geodesies defined as energy minimizing paths gives rise to a discrete logarithmic map, a variational definition of a 
discrete exponential map, and a time discrete parallel transport. This new concept is applied to a shape space in which 
shapes are considered as boundary contours of physical objects consisting of viscous material. The flexibility and 
computational efficiency of the approach is demonstrated for topology preserving shape morphing, the representation 
of paths in shape space via local shape variations as path generators, shape extrapolation via discrete geodesic flow, 
and the transfer of geometric features. 
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1. Introduction. Geodesic paths in shape space allow to define smooth and in some 
sense geometrically or physically natural connecting paths 0{t), t G [0, 1], between two 
given shapes 0(0) , or they enable the extrapolation of a path from an initial shape 0(0) 
and an initial shape variation SO which encodes the path direction. Applications include 
shape modeling in computer vision [17, 16], computational anatomy, where the morphing 
path establishes correspondences between a patient and a template [2, 26], shape clustering 
based on Riemannian distances [32], as well as shape statistics [9, 13], where geodesic paths 
in shape space transport information from the observed shapes into a common reference frame 
in which statistics can be performed. 

As locally length minimizing paths, geodesic paths require to endow the space of shapes 
with a Riemannian metric which encodes the preferred shape variations. There is a rich 
diversity of Riemannian shape spaces in the literature. Kilian et al. compute isometry invari- 
ant geodesies between consistently triangulated surfaces [16], where the Riemannian metric 
measures stretching of triangle edges, while the metric by Liu et al. also takes into account 
directional changes of edges [22]. 

For planar curves, different Riemannian metrics have been devised, including the L^- 
metric on direction and curvature functions [17], the -metric on stretching and bending 
variations [31], as well as curvature-weighted L^- or Sobolev-type metrics [25, 34], some of 
which allow closed-form geodesies [37, 33]. A variational approach to the computation of 
geodesies in the space of planar Jordan curves has been proposed by Schmidt et al. in [30]. 
The extrapolation of geodesies in the space of curves incorporating translational, rotational, 
and scale invariance has been investigated by Mennucci et al. [24]. A Riemannian space of 
non-planar elastic curves has very recently been proposed by Srivastava et al. [32]. 

In the above approaches, the shape space is often identified with a so-called pre-shape 
space of curve parameterizations over a ID domain (or special representations thereof) mod- 
ulo the action of the reparameterization group. It is essential that the metric on the pre-shape 
space is invariant under reparameterization or equivalently that reparameterization represents 
an isometry in the pre-shape space so that the Riemannian metric can be inherited by the 
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Fig. 1.1. Top: Given the first shape on the left and an initial variation ("i described by the second shape, 
a discrete geodesic path is extrapolated. Bottom: The texture of a video frame can be transported along with the 
resulting geodesic flow. 

shape space. Such reparameterization-invariant metrics can also be defined on the space of 
parameterized 2D surfaces [1, 20]. For certain representations of the parameterization one is 
lead to a very simple form of the metric, e.g. an L^-type metric [19]. 

The issue of reparameterization invariance does not occur when the mathematical de- 
scription of the shape space is not based on parameterizations, which often simplifies the 
analysis (and is also the approach taken here). When warping objects in R'^, a shape tube 
in IR^+^ is formed. Zolesio investigates geodesic in terms of shortest shape tubes [39]. The 
space of sufficiently smooth domains (9 C can be assigned a Riemannian metric by iden- 
tifying the tangent space at O with velocity fields v : O and defining a metric on these. 
Dupuis et al. employ a metric 

Q{v.,v) = I Lv ■ V dx 
Jd 

for a higher order elliptic operator L on some computational domain D C [7], ensuring a 
diffeomorphism property of geodesic paths. A corresponding geodesic shooting method has 
been implemented in [3]. Fuchs et al. propose a viscous-fluid based Riemannian metric [12]. 
Fletcher and Whitaker employ a similar metric on pullbacks of velocity fields onto a reference 
shape [10]. Miller and Younes consider the space of registered images as the product space of 
the Lie group of diffeomorphisms and image maps. They define a Riemannian metric using 
sufficiently regular elliptic operators on the diffeomorphism-generating velocity fields, which 
may also depend on the current image [27] . A morphing approach based on the concept of 
optimal mass transport has been proposed by Haker et al. [14, 38]. An image or a shape is 
viewed as mass density, and for two such densities po, Pi : ^ ^ R the Monge-Kantorovich 
functional 

\ip{x) — x\'^po{x) dx 

is minimized over all mass preserving mappings i/j : D ^ D, i.e. mappings with po = 
pi o ijjdeiVijj. A morphing path then is given by p{t) = po o det VV^(t)~^ for 

^(t) = tip -\- {1 — t)id, t G [0, 1]. Like for our approach there is a continuum-mechanical 
interpretation of minimizing the action of an incompressible fluid flow [4], however, the flow 
typically does neither preserve local shape features or isometrics nor the shape topology. 

Very often, geodesies in shape space are approached via the underlying geodesic evolu- 
tion equation, and geodesies between two shapes are computed by solving this ODE within a 
shooting method framework [17, 3, 1]. An alternative approach exploits the energy-minimi- 
zing property of geodesies: Schmidt et al. perform a GauB-Seidel type fixed-point iteration 
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which can be interpreted as a gradient descent on the path energy, and Srivastava et al. de- 
rive the equations of a gradient flow for the path energy which they then discretize [32]. In 
contrast, we employ an inherently variational formulation where geodesies are defined as 
minimizers of a time discrete path energy. Discrete geodesies are then defined consistently as 
minimizers of a corresponding discrete energy. 

In this paper we start from this time discretization and consistently develop a time dis- 
crete geodesic calculus in shape space. The resulting variational discretization of the basic 
Riemannian calculus consists of an exponential map, a logarithmic map, parallel transport, 
and finally an underlying discrete connection. To this end, we replace the exact, computation- 
ally expensive Riemannian distance by a relatively cheap but consistent dissimilarity measure. 
Our choice of the dissimilarity measure not only ensures consistency for vanishing time step 
size but also a good representation of shape space geometry already for coarse time steps. 
For example, rigid body motion invariance is naturally incorporated in this approach. We 
illustrate this approach on a shape space consisting of homeomorphic viscous-fluid objects 
and a corresponding deformation-based dissimilarity measure. 

Different from most approaches, which first discretize in space via the choice of a param- 
eterization, a set of control points, or a mesh, and then solve the resulting transport equations 
by suitable solvers for ordinary differential equations (see the discussion above), our time 
discretization is defined on the usually infinite dimensional shape space. It results from a 
consistent transfer of time continuous to time discrete variational principles. Thereby, it leads 
to a collection of variational problems on the shape space, which in our concrete implemen- 
tation of the proposed calculus consists of non-parameterized, volumetric objects. 

Let us also already mention a further remarkable conceptual difference. The way the time 
discrete geodesic calculus is introduced differs substantially from the way the time continuous 
counterpart is usually developed. In classical Riemannian differential geometry one first 
defines a connection (v^w) VyW for two vector fields v and w on a. manifold Ai. With 
the connection at hand a tangent vector w can be transported parallel along a path with motion 
field V solving = 0. Studying those paths where the motion field itself is transported 
parallel along the path (Le. it solves the ODE VyV = 0) one is led to geodesies. Next, the 
exponential map is introduced via the solution of the above ODE for varying initial velocity. 
Finally, the logarithm is obtained as the (local) inverse of the exponential map. 

In the time discrete calculus we start with a time discrete formulation of path length and 
energy and then define discrete geodesies as minimizers of the discrete energy. Evaluating 
the initial step of a discrete geodesic path as the discrete counterpart of the initial velocity we 
are led to the discrete logarithm. Then, the discrete exponential map is defined as the inverse 
of the discrete logarithm. Next, discrete logarithm and discrete exponential allow to define a 
discrete parallel transport based on the construction of a sequence of approximate Riemannian 
parallelograms (commonly known as Schild's ladder [8]). Finally, with the discrete parallel 
transport at hand, a time discrete connection can be defined. 

Let us note that the approximation of parallel transport in shape space via Schild's ladder 
has also been used in the context of the earlier mentioned flow of diffeomorphism approach 
[29, 23]. In our discrete framework, however, the notion of discrete parallel transport is 
directly derived from the parallelogram construction, consistently with the overall discrete 
approach to geodesies. 

A related approach for time discrete geodesies has been presented in an earlier paper by 
Wirth et al. [36]. In contrast to [36], we here do not restrict ourselves to the computation 
of geodesic paths between two shapes but devise a full-fledged theory of discrete geodesic 
calculus (cf. Figure 1). Furthermore, different from that approach we ensure topological con- 
sistency and describe shapes solely via deformations of reference objects instead of treating 
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deformations and level set representations of shapes simultaneously as degrees of freedom, 
which in turn strongly simplifies the minimization procedure. 

The paper is organized as follows. In Section 2 we introduce a special model for a shape 
space, the space of viscous fluidic objects, to which we restrict our exposition of the geodesic 
calculus. Here, in the light of the discrete shape calculus to be developed, we will review the 
notion of discrete path length and discrete path energy. After these preliminaries the actual 
time discrete calculus consisting of a discrete logarithm, a discrete exponential and a discrete 
parallel transport together with a discrete connection is introduced and discussed in Section 
3. Then, Section 4 is devoted to the numerical discretization via characteristic functions and 
a parameterization via deformations over reference paths. Finally, we draw conclusions in 
Section 5. 

2. A space of volumetric objects and an elastic dissimilarity measure. To keep the 
exposition focused we restrict ourselves to a specific shape model, where shapes are repre- 
sented by volumetric objects which behave physically like viscous fluids. In fact, the scope 
of the variational discrete geodesic calculus extends beyond this concrete shape model. We 
refer to Section 5 for remarks on the application to more general shape spaces. 

2.1. The space of viscous-fluid objects. Let us introduce the space M. of shapes as the 
set of all objects O which are closed subsets of (d = 2, 3) and homeomorphic to a given 
regular reference object O^ef, i-e. O = ^(Oref) for an orientation preserving homeomorphism 
(j). Furthermore, objects which coincide up to a rigid body motion are identified with each 
other. A smooth path (0(t))t^[o,i] in this shape space is associated with a smooth family 
(0(t))^^[o,i] of deformations. To measure the distance between two objects, a Riemannian 
metric is defined on variations 50 of objects O ^ M which reflects the internal fluid friction 
— called dissipation — that occurs during the shape variation. The local temporal rate of 
dissipation in a fluid depends on the symmetric part e[v\ := ^ {Vv + Vv^) of the gradient of 
the fluid velocity v : O ^ IR^ (the antisynmietric remainder reflects infinitesimal rotations), 
and for an isotropic Newtonian fluid, we obtain the local rate of dissipation 

diss(V'?;) = \{tYe[v\f + 2/itr(eH^) , (2.1) 

where A, /i are material- specific parameters. Given a family ((/)(t))^^[o,i] of deformations of 
the reference object O^ef. the change of shape along the path (O(t))^^[o,i] can be described by 
the (Lagrangian) temporal variation ^{t) or the associated (Eulerian) velocity field 

v{t) = m o <t>-\t) 

on O. Hence, the tangent space TqM to Ai 3.1 a. shape O can be identified with the space 
of initial velocities v = 0(0) o 0~^(O) for deformation paths with 0(0, O^ef) = O. Here we 
identify those velocities v which lead to the same effective shape variation, i.e. those with the 
same normal component v • non dO, where n is the outer normal on dO. Now, integrating 
the local rate of dissipation for velocity fields v onO = 0(0, O^ef), we define the Riemannian 
metric Go on TqM. as the symmetric quadratic form with 

Go{v,v)= min / diss(Vi;(x)) dx . (2.2) 

{v \ v-n=v-n on dO} Jq 

For the shape variation along a path O : [0, 1] ^ described by the Eulerian motion field 
(^(t))t^[o,i], path length L and energy E are defined as 

L[(OW)te[o,i]] = /o ^QoitMt),v{t))^t , (2.3) 
E[(O(i))te[0,i]] = jlQoitMt),v{t))At. (2.4) 
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Fig. 2.1. Discrete geodesic between the letters Q and A, extrapolated beyond A. Colors indicate the local rate 
of dissipation (from blue, low, to red, high). 

Paths which (locally) minimize the energy E or equivalently the length L are called geodesies 
(c/. Figure 2.1). A geodesic thus mimics the energetically optimal way to continuously de- 
form a fluid volume. 

2.2. Approximating the distance. The evaluation of the geodesic distance based on a 
direct space and time discretization of (2.2) and (2.3) turns out to be computationally very 
demanding (c/. for instance the approaches in [3, 7]). Hence, we use here an efficient and 
robust time discrete approximation based on an energy functional W which locally behaves 
like the squared Riemannian distance {i.e. the squared length of a connecting geodesic): 

Given two shapes O and (5, we consider an approximation 

dist'((9,(5) ^ WoM, (2.5) 

where Wo [^] '= Jq dx is the stored deformation energy of a deformation : O ^ 

and ip is the minimizer of this energy over all such deformations with = O. Here, 

: R^'^ ^ [0, oo) is a so-called hyperelastic energy density. 

In correspondence to our assumption that objects are identical if they coincide up to a 
rigid body motion, we require W to be rigid body motion invariant. Furthermore, we assume 
the objects to have no preferred material directions so that W is in addition isotropic, which 
altogether leads to W{RAU) = W{A) for aU R,U e SO{d),A e R^'^ (cf. [5]). In 
the undeformed configuration for i/j = id, energy and stresses (the first derivatives of W) 
Sire supposed to vanish so that we require W{1) = 0, VW{1) = (where VW denotes 
the derivative with respect to the matrix argument). Furthermore, we need oo 
as det ^4 ^ to prohibit material self-penetration, which is linked to the preservation of 
topology. The approximation property (2.5) relies on a consistent choice of Wo for the given 
metric Qo which can be expressed by the relation 



= 11 V^W{l){\/v,\/v)dx= [ diss(Vv)dx (2.6) 
=0 ^ Jo Jo 



along any object path 0{t) = ilj{t,0),t G R, with ?/;(0) = id and velocity field v = ^(0). 
Using the notion of the Hessian Hess^i of a function on a manifold as the endomorphism 
representing its second variation in the metric, we can rephrase this approximation condition 
more geometrically as 

^Hess^W(9[id] id 

with the usual identification of objects O and deformations ^. For the deformation energy 
density W, this condition implies that its Hessian V'^W{1) has to satisfy ^V'^W{1){A, A) = 
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diss(A) for all A e M.^'^ . A suitable example is 



l^(A) = ^tr(AM) + -detA^- (^/i+-JlogdetA ^ ^. 

Assume that the energy density satisfies the above-mentioned properties. We observe that the 
metric Qo is the first non- vanishing term in the Taylor expansion of the squared length of a 
curve, i.e. 

{h[{0{t)\^^,,T]]f =T^Qoi0)(.v,v)+O{T') 

with V = ^(0) o ^~^(0) being the initial tangent vector along a smooth path (0(t))^^[o,T] = 
(^(t, (9ref))tG[o,i] • Thus, since the Hessian of the energy Wo and the metric Qq are related by 
(2.6), we obtain that the second order Taylor expansions of dist^(0, ^(O)) and >Vc)[^] in 
coincide and indeed 

dist^{0, 0) = min Wo[^] + 0(dist^(C>, O)) . (2.7) 

Here, different from [36] we neither take into account mismatch penalties nor perimeter reg- 
ularizing functionals for each object O/^, = 0, i^. 

2.3. Discrete length and discrete energy. Now, we are in a position to discretize length 
and energy of paths ((9(t))te[o,i] in shape space. To this end, we first sample the path at times 
tk = kr for k = 0, K (r = ^), denote Ok := 0{tk), and obtain the estimates 

mO{t))te[o,i]] > Ef=idist(Ofe_i,Ofc) 

E[((!?W)te[o,i]] > ^Ef=idist'(Ofc-i,Ofe) 

for the length and the energy, where equality holds for geodesic paths. Indeed, the first 
estimate is straightforward, and the application of the Cauchy-Schwarz inequality leads to 



^dist2(Ofc_i,Ofc) <J2\ [ \/Go(tMt).v{t))dt] 

K r^kr 

<E^ / ^o(t)(^W,^W)clt = TE[(O(t)),e[0,i]] 



which implies the second estimate. 

Together with (2.7) this motivates the following definition of a discrete path energy and 
a discrete path length for a discrete path (Oq? • • • ? Ok) in shape space: 

L[((9o, . . . , Ok)] = Ef=i V^o,_A^k] , (2.8) 
E[((9o, . . . , Ok)] = \ Ek=i m , (2.9) 

where V^/, = argmin|^ | ^i^Ok-i)=Ok} "^o^-i M (^/- also [36]). In fact, (2.8) and (2.9) can 
for general smooth paths even be proven to be first order consistent with the continuous 
length (2.3) and energy (2.4) as r ^ 0. For illustration, if is a two-dimensional manifold 
embedded in R^, we can interpret the terms Wo^.i as the stored elastic energies in springs 
which connect a sequence of points Ok on the manifold through the ambient space. Then the 
discrete path energy is the total stored elastic energy in this chain of springs. 
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Fig. 2.2. Nonlinear video interpolation via a discrete geodesic (top) between two segmented photographs of 
white and red blood cells (first and last picture of bottom row, courtesy Robert A. Freitas, Institute for Molecular 
Manufacturing, California, USA). The bottom row shows pushforwards and pullbacks of the end images under the 
deformations along the discrete geodesic. 




A discrete geodesic (of order K) is now defined as a minimizer of E[((9o, • • • , Ok)] for 
fixed end points (9o, Ok- The discrete geodesic is thus an energetically optimal sequence of 
deformations from Oq into Ok- 

In the minimization algorithm to be discussed in Section 4.1 we do not explicitly mini- 
mize E[(Oo, • • • , Ok)] for the object contours as in [36] but instead for reference deforma- 
tions defined on fixed reference objects. Figure 2.2 shows a discrete geodesic in the context 
of multicomponent objects, which is visually identical to that obtained by the more com- 
plex approach in [36]. Here, deformations are considered which map every component of a 
shape onto the corresponding component of the next shape in the discrete path as the obvious 
generalization of discrete geodesies between single component shapes. 

While in the continuous case geodesic curves equally min- 
imize length (2.3) and energy (2.4), minimizers of the discrete 
path length (2.9) are in general not related to discrete geodesies 
(and thus also not to continuous geodesies as r ^ 0). Indeed, 
let us consider a two-dimensional manifold A4 embedded in 
R^, paired with the deformation energy W^fc-Jid + Cfe] •= 

IC/cP for a displacement vector Cfc in connecting points Fig. 2.3. A continuous 
Ok-i and Ok on M. Now take into account a continuous geodesic and a discrete path which 

almost minimizes the discrete path 

geodesic and a discrete path on M where the end points are /^^^^/, ^„ ^ two-dimensional man- 
close to each other in the embedding space but far apart on if old embedded in . 
the surface. Figure 2.3 depicts such a configuration with a dis- 
crete path which almost minimizes the discrete path length. A minimizer of the discrete 
path length will always jump through the protrusion and never approximate the continuous 
geodesic, whereas minimizers of the discrete path energy satisfy Wc)^_Jid + Cfc] ^0 as 
r ^ and thus rule out such a short cut through the ambient space. 

3. Time discrete geodesic calculus. With the notion of discrete geodesies at hand we 
will now derive a full-fledged discrete geodesic calculus based on a time discrete geometric 
logarithm and a time discrete exponential map, which then also give rise to a discrete parallel 
transport and a discrete Levi-Civita connection on shape space. 

3.1. Discrete logarithm and shape variations. If ((9(t))^^[o,i] is the unique geodesic 
on M connecting (9 = (9(0) and O = 0(1), the logarithm of O with respect to O is defined 
as the initial velocity v ^TqM of the geodesic path. In terms of Section 2.1 we have 

\ogo{d) = v{()) 

for v{t) = (j){t) o (j)(t)~^, where ^(t, 0,^f) = 0{t) defines the associated family of defor- 
mations. On a geodesically complete Riemannian manifold the logarithm exists as long as 
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dist((9, O) is sufficiently small. The associated logarithmic map log^-, : O ^ v{0) G TqM 
represents (nonlinear) variations on the manifold as (linear) tangent vectors. 

The initial velocity ^'(0) can be approximated by a difference quotient in time, 

^(O,x) = -C(x) + O(r), 

T 

where ({x) = (/)(r, x) o 0(0, x)~^ — x denotes a displacement on the initial object O. Thus, 
we obtain 

rlogc(a) = C(x) + 0(r2). 

This gives rise to a consistent definition of a time discrete logarithm. Let (Oq, . . . , Ok) be a 
discrete geodesic between O = Oq and O = Ok with an associated sequence of matching 
deformations . . . , i/jk, then we consider ^^i for the displacement Ci(^) = V^i(^) — ^ as 
an approximation of ^'(O) = logQ{0). Taking into account that r = we thus define the 
discrete -logarithm 

(^LOG)^(a) := Ci . (3.1) 
In the special case K = 1 and a discrete geodesic ((9, O) we simply obtain 
{\LOG)^{d) = argmin^^^ , (id+Ci)(0)=a} ^o[^d + Ci] • 

As in the continuous case the discrete logarithm can be considered as a representation of 
the nonlinear variation O of O in the (linear) tangent space of displacements on O. On a 
sequence of successively refined discrete geodesies we expect 

K{^LOG)^{d) ^ logc(a) (3.2) 

for K ^ (X) (cf. Figure 3.1 for an experimental validation of this convergence behaviour). 

3.2. Discrete exponential and shape extrapolation. In the continuous setting, the ex- 
ponential map exp^ maps tangent vectors v e TqM onto the end point 0{1) of a geodesic 
((9(t))t^[o,i] with (9(0) = O and the given tangent vector v at time 0. Using the notation 
from the previous section we have expQ{v{0)) = O and, via a simple scaling argument, 
expQ = 0{^) for k = 0, . . . , K. We now aim at defining a discrete power k 

exponential map EXPq such that EXP|,(Ci) = Ok on a. discrete geodesic {O, Oi, . . . , Ok) 
of order K > k with ("i = {^LOG) ^{O) (the notation is motivated by the observation that 
exp(A:s) = exp^(s) on R or more general matrix groups). Our definition will reflect the 
following recursive properties of the continuous exponential map, 

expo(l^) = {\\ogo)~^ (v), 
expo(2v) = (^log^)"^ {v), 
expo(to) =exp^^p^((;,_2)^)(2v/e-i) 

for Vk-i := logexp^((/c-2)^) expc)((A: - 1)^) . 
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Fig. 3.1. In the bottom left we experimentally verify the convergence stated in (3.9) by computing discrete 
geodesies starting from 'P' and ending at EXPp{C/K) for K = 1,2, 3, 4, 8 and C = 8(|L0G)^(A). In the 
top right, based on a computation of discrete geodesies between 'P' and 'A' of order K = 8,4, 3, 2, 1, the resulting 
logarithms (-^LOG)^(A) are depicted, where (-^LOG)^(A) is the displacement from 'P' to the second shape 
of the respective geodesic. To the left of the initial shape 'P', those shapes are displayed which result from applying 
the displacement K(^LOG) ^(A) to 'P' to experimentally verify convergence as stated in (3.2). The arbitrary 
rotations are due to the rigid body motion invariance of our approach. 



Replacing exp(A:-) by EXP^, \ log by (^LOG), and \ log by (^LOG) we obtain the recur- 
sive definition 



EXP^(C) 
EXP2,(C) 
EXP^(C) 



=(iLOG)-'(C), (3.3) 
= (iLOG)-'(C), (3.4) 

withCfe-i := (iLOG)gxp^_.(^)EXP^-i(C). 



It is straightforward to verify that EXP^ = (^LOG)^^ as long as the discrete logarithm 
on the right is invertible. Equation (3.3) implies EXP^(C) = (id + 0{O), and (3.4) in fact 
represents a variational constraint for a discrete geodesic flow of order 2 : 

Given the object O we consider discrete geodesic paths (O, Oi, ^^2) of order 2, where 
for any chosen O2 the middle object Oi is defined via minimization of (2.9) so that we may 
write di[d2\ We now identify EXP as the object 62 for which (id + C)(0) = Oi[^2], 
i.e. C is the energetically optimal displacement from O to Oi [O2] and thus satisfies 

id + C = argmin^^^ |^^(^)^^^[^^]-^ Woi^Pi] (3.6) 
up to a rigid body motion. 

Alternatively, the condition (3.6) can be phrased as 

id + C = argmin^^^} ^^{^^2 1 (^20^iKO)=62} {^o[^i] + y\^MO)M • (3.7) 
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Fig. 3.3. Left: Discrete geodesic between two end shapes, where color indicates the local rate of dissipation 
(from blue, low, to red, high). Right: Given the first shape and its variation (in terms of the optimal matching 
deformation to the second shape), a discrete geodesic is extrapolated. 



Figure 3.2 conceptually sketches the 
procedure to compute EXP|,(C). For given 
initial object O and initial displacement C 
the discrete exponential EXP|,(C) is se- 
lected from a fan of discrete geodesies with 
varying O2 as the terminal point of a dis- 
crete geodesic of order 2 in such a way that 
(3.7) holds. To compute EXP|,(C) in the 
geodesic flow algorithm (3.4) and (3.5) we 
have to find the root of 




Fig. 3.2. Conceptual sketch of the procedure to 
compute EXP^ {C,) with objects represented as points. 



i^o,c(^2) = (^LOG)^((92)-C, (3.8) 

implicitly assuming that C is small enough so that discrete geodesies are unique (c/. Sec- 
tion 4.2 for the algorithmic realization based on a representation of the unknown domain O2 
via a deformation). 

Equation (3.5) describes the recursion to compute EXP^ based on the above EXP^ sin- 
gle step scheme: Forgiven Ok-2 = EXP^~^(C) and Ok-i = EXP^~^(C) one first retrieves 
= — id from the previous step, where 

tljk-i = argmini^ , ^(o,_2)=o,_i} M • 

Then (3.4) is applied to compute Ok from Ok-2 and (k-i as the root of Fof^_2,Ck-i- 

For sufficiently small ^ we expect EXP^ to be well-defined. Since by definition, every 
triplet ((9/e_i,(9/c,(9/c+i)ofthe sequence Ok = EXPq{() is a geodesic of order 2 and mini- 
mizes E[((9/c_i, O/c, Ok-\-i)], the resulting family {Ok)k=o,...,K indeed is a discrete geodesic 
of order K. In fact, discrete geodesies that are variationally described as discrete energy min- 
imizing paths between two given objects can be reproduced via the discrete geodesic flow 
associated with the discrete exponential map (cf. Figure 3.3). 

As for the discrete logarithm we experimentally observe convergence of the discrete 
exponential map in the sense 

EXP^ (K) ^ exp^(C) for ^ 00 (3.9) 

as shown in Figure 3.1. An example of geodesic shape extrapolation for multicomponent 
objects is depicted in Figure 3.4. 

3.3. Discrete parallel transport and detail transfer. Parallel transport allows to trans- 
late a vector ( e ToM. (which is considered as the variation of an object O = O{0)) along 
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Fig. 3.4. Top: Given the first shape and its initial variation Ci represented by the second shape, a discrete 
'esic is extrapolated. Bottom: The texture of a video frame can be transported along with the varying shapes. 




a curve (O(t))t^[o,i] in shape space. The resulting (C(t))t^[o,i] changes as little as possible 
while keeping the angle between C,{t) and the path velocity v{t) fixed. Using the Levi-Civita 
connection this can be phrased as V^(t)C(^) = 0. There is a well-known first-order approxi- 
mation of parallel transport called Schild's ladder [8, 15], which is based on the construction 
of a sequence of geodesic parallelograms, sketched in Figure 3.5, where the two diagonal 
geodesies always meet at their midpoints. Given a curve (O(t))^^[o,i] and a tangent vector 
C/c-i G T'(9((/p_i)^)A^, the approximation C,k ^ ^c(/cr)-^ of the parallel transported vector 
via a geodesic parallelogram can be expressed as 

Here, is the midpoint of the two diagonals of the geodesic parallogramm with vertices 
0{{k — l)r), Ol_^, O^., and 0{kr). This scheme can be easily transferred to discrete curves 
((9o, . . . , Ok) in shape space based on the discrete logarithm and the discrete exponential 
introduced above. In the kth step of the discrete transport we start with a displacement (k-i 
onOk-i and compute 
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Fig. 3.6. Discrete parallel transport is applied along two discrete geodesic paths connecting letters 'P' and 
'A' (left) and two different poses of a dog (right), respectively. On the top the original discrete geodesies are shown, 
while the bottom left shows the discrete parallel transport of serifs along the geodesic between the two letters, and 
the bottom right shows the transport of changes on the first dog 's shape, which allows to copy the changes to the 
other poses. 



k-1 


— EXPo^_^Cfe-i > 






= EXP^P_^ ((^LOG)^, 




Ol 


= EXP^._. ((tLOG)^^ 


jo^)) 




= (tLOG)^joD, 





where (k is the transported displacement on Ok • Here, is the midpoint of the two discrete 
geodesies of order 2 with end points 0^_^, Ok and Ok-i, O^., respectively. Since the last of 
the above steps is the inverse of the first step in the subsequent iteration, these steps need to 
be performed only for k = K. We will denote the resulting transport operator by Pok,---,Oo • 
Figure 3.6 shows examples of discrete parallel transport for feature transfer along curves in 
shape space. 

Remark: As in the continuous case, the discrete parallel transport can be used to define 
a discrete Levi-Civita connection. For ^ eToM and for a vector field r] in the tangent bundle 
TA4 one computes Or = EXP^(r^) and then defines 

V^rj := - {Po,oMOr) - viO)) 

S J. 

as the time discrete connection with time step size r. 

4. Numerical discretization. The proposed discrete geodesic calculus requires an ef- 
fective and efficient spatial discretization of 

- volumetric objects O in the underlying shape space, 

- of nonlinear deformations to encode matching correspondences, 

- and of linear displacements ( as approximate tangent vectors. 

We restrict ourselves here to the case of objects (9 C R^. To this end we consider the space 
Vh of piecewise affine finite element functions on a regular simplicial mesh over a rectangular 
computational domain D. Here h indicates the grid size, where in our applications h ranges 
from a coarse grid size to a fine grid size 2~^. Then, deformations and displacements are 
considered as functions in (V/i)^. Objects O, the original degrees of freedom in our geometric 
calculus, will be represented via deformations (j) over reference objects O (e.g. 0,^f), i.e. 
O = (j){0). These reference objects are encoded by approximate characteristic functions 
G Vh and the deformations (j) are considered as injective deformations (j) : D ^ and 
discretized as elements in (V/i)^. 

4.1. Parameterization of discrete geodesies. To compute a discrete geodesic — dif- 
ferent from [36] — we now replace the objects (9o, . . . , Ok as arguments of the energy (2.9) 
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Fig. 4.1. A diagram illustrating the parameterization of the domains along the discrete geodesic via 
deformations (f)k — indicated in red as the actual degrees of freedom — over reference domains Ok- 



by associated deformations 0o , . . . , 0k over a set of reference domains Oq , • • • , Ok as de- 
scribed above. By this technique, instead of K deformations and K — 1 domain descriptions 
(e.g. via level sets) as in [36] we will be able to consider solely K -\- 1 parameterizing de- 
formations, which turns out to be a significant computational advantage. Next, we assume 
that reference matching deformations . . . , iIjk are given with Ok = tpk{Ok-i) (cf. Fig- 
ure 4.1). Now, we express the matching deformations {ipk)k=i,...,K over which we minimize 
in (2.9) in terms of the parameterizing deformations (0/c)fe=o,...,K and the reference matching 
deformations {tl^k)k=i,...,K and set 

i^k = (t>k ^ i^k ^ (t>k-i 

for = 1, i^r. Now, using a change of variables one can rewrite the deformation energies 
Wo,_AM in (2.9) as 

X6,_W{V{^k o ^fc)(V0fc_i)-^) det V0fc_idf . 



D 



Furthermore, instead of void we consider a 5i times softer material outside the object domain 
replacing Xq = (1 — ^^i)Xc) +^i so that altogether the deformation energy Wok-i [V^/c] 

is replaced by the following energy over the parameterizing deformations 0^-1 and 0/^: 



W^-^'^[0fe_i,0fe]= / x'^ iy(V(0feO^,)(V0,_i)-^)detV0,_idx 



The condition that Oq and Ok are prescribed is taken care of with penalty functional 

^6%- [^^] " " / ^^^^ * ^Oi ~ * XOi o 0i)^ 

for i = O^K, where Gs^ is a Gaussian of filter width 62. Finally, to ensure that not only the 
concatenations 0/c o ?/;^ o 0^1 ^ of deformations are regular but also every single deformation 
0/c, we add a term ^3Wd[0/c] for all /c = 0, . . . , iC. Summarizing, to compute a discrete 
geodesic between two shapes Oq , Ok or a discrete logarithm, we minimize the total energy 

fe=l k=0 

over all the parameterizing deformations 0o, . . . , 0k . This minimization then determines the 
shapes Ok = (})k{Ok) forming the discrete geodesic as well as the discrete logarithm 

{j^LOG) {Ok) = ^1 0^1 o - id. 
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Fig. 4.2. Discrete geodesies with three, five, and nine shapes (top to bottom). Despite large shape variations 
and strong local rotations, the choice of the nonlinear energy W yields good results already for coarse time steps. 

Due to the assumptions on the energy integrand W, self-penetration of deformed objects is 
ruled out in our fully discrete model, which as a byproduct leads to topology preservation 
along discrete geodesies. The resulting matching deformations ^/^i, . . . , ijjK are determined 
up to translations and rotations due to the built-in frame indifference (c/. Section 2.2). In 
our applications the parameters of the algorithm are chosen as follows: (5i = 0.01, (^2 = /i, 

^3 =0.01,6 = 0.1. 

Furthermore, the previously described finite element discretization is applied, and the 
energies are computed via Simpson quadrature on each element. As opposed to Gaussian 
quadrature, this quadrature rule has the advantage that quadrature points lie in the corners 
of each finite element, which is where the extremal values of the finite element functions 
and their gradients occur. This is relatively important when dealing with deformations as 
finite element functions, since a Gaussian quadrature rule might for example miss that a 
deformation exhibits self-penetration in the corner of one element. Such deformations will 
thus not be rejected during energy minimization, which in turn results in technical difficulties 
when computing pullbacks or pushforwards with respect to these deformations and when 
prolongating them onto grids with finer resolution. 

Pullbacks of a finite element function f GVh with respect to a finite element deformation 
(j) G (V/i)^ at quadrature points x are computed by first evaluating the deformation at that 
point. lf(j){x) lies outside the computational domain I), it is projected back onto the boundary 
dD (this approach is more robust than e.g. neglecting such points, since it avoids structural 
changes when deformed quadrature points toggle between inside and outside). Finally, / is 
evaluated at the projected position. 

The minimization is performed by a Newton trust region algorithm [6] on the (i^T + 1)- 
tuple of deformations (^o, • • • , <I^k) and requires the evaluation of first and second derivatives 
of the energy. Note that the second derivative involves mixed derivatives with respect to 
(j)k-i and for all /c = 1, i^T. During the Newton iteration, these mixed terms provide 
a coupling between all deformations which results in a fast relaxation and balance between 
them. 

To improve the efficiency of the resulting method, a cascadic finite element approach 
is used which proceeds from a coarse to a fine resolution of objects and deformations on a 
dyadic hierarchy of meshes. Simultaneously the resolution of the discrete geodesies can be 
increased in time (cf. Figure 4.2). In case very large nonlinear deformations occur during 
the optimization, from time to time the reference objects (Oq, • • • , Ok) are replaced by the 
current object approximations ((9o, . . • , Ok) and the deformations (j)k are reset to the identity 
deformation id for /c = 0, i^T. 

4.2. Solving the constrained optimization problem associated with EXP^. To find 
the zero of (3.8) we employ the optimality condition associated with the variational definition 
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of (|L0G). Indeed, for given ( we introduce two deformations ^i, ^2 and corresponding 
objects O2 := {ip2 o (id + C)) ^1 •= V^i(C^)- Now, we associate to the discrete path 
(0, with underlying deformations ipijip2 the energy 

£[^1,^2] [ W{\/^i)dx^ [ l^(V(7/;2o(id + C)o^r')) dx. 
Jo Jipi(0) 

We obtain S^^f [V^i , 1^2] = as the necessary condition to ensure that Sl^pi , ?/^2] actually rep- 
resents the path energy in (2.9) connecting the two objects O and O2 via a discrete geodesic 
of order 2 with the intermediate object ipi{0). In the notion of Section 3.2 the property of 
(O, (5i , 02) to be a discrete geodesic can be phrased as ipi (O) = Oi [O2]. Thus, a necessary 
condition for (3.7) to hold is given by the condition 

a^.f [7/^1,7/^2] Ui=id+c = (4.1) 

for the remaining unknown 1JJ2. With respect to the algorithmic realization we reformulate 
and regularized the energy as described in Section 4.1 above to yield 

£^'[^i,H = I Xo' (^(V^i) + W^(V(i/^2o(id + C))(V^i)-')detVi/^i) dx. 

JD 

Now, the discrete counterpart of (4.1) is the condition 

^=d^,£'^[^^,^2]\^^^,^^^ {0) 

- VW(V^2 o (id + C)) : iy^2 o (id + Q) VO (1 + VC)"' det(l + VQ 
+ W{yiJ2 o (id + C)) det(l + VC) tr((l + VC)"^ V6>)) dx , (4.2) 

where A : B = ti^A^ B) for matrices B G R^'^. This equation has to hold for all test 
deformations 0. In our finite element context, the corresponding test functions are taken to be 
all finite element basis functions so that (4.2) becomes a system of nonlinear equations which 
is solved for ^^2 via Newton's method. Here too, we first find ip2 on a coarse grid and then 
use the resuk as initiahzation of Newton's method on finer grids. 

5. Conclusions and outlook. Based on a variational time discretization of geodesic 
paths in shape space we have proposed a novel time discrete geodesic calculus, which consists 
of discrete logarithmic and exponential maps, discrete parallel transport and a discrete con- 
nection. We demonstrate how to use this discrete calculus as a robust and efficient toolbox for 
shape morphing, shape extrapolation, and transport of shape features along paths of shapes. 
Although in this expository article we restricted ourselves to two-dimensional objects, the 
approach can be carried over to 3D viscous-fluid shapes. The concept can also be adapted to 
deformations of hypersurfaces and corresponding deformation energies, which measure tan- 
gential as well as normal bending stresses [11,21]. For example, a generalization to the space 
of planar elastic curves and thin shell surfaces is feaible. Furthermore, instead of a metric 
structure induced by the viscous flow paradigm, the Wasserstein distance of optimal transport 
can be considered [35]. It this case the time discretization of the Monge-Kantorovich prob- 
lem proposed by Benamou and Brenier [4] is a possible starting point. Beyond these future 
directions of generalization, a theoretical foundation has to be established with existence and 
regularity results for the above-mentioned infinite dimensional shape spaces. Furthermore, 
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the limit behaviour of the discrete geodesic calculus for vanishing time step size and the con- 
vergence of the discrete path energy to the corresponding continuous path energy in the sense 
of F-convergence has to be investigated (c/. the work by Miiller and Ortiz on F-convergence 
of a time discrete action functional in the case of Hamiltonian systems [28]). Finally, given 
the notion of a time discrete transport, the relation of the curvature tensor to the parallel trans- 
port along the edges of a quadrilateral (c/. Proposition 1.5.8. in [18]) can be used to define 
a time discrete curvature tensor, which then allows an exploration of the local geometry of 
shape space. 
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